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We present a new model of surface transport in premelted films that is applicable to a wide 
range of materials close to their melting point. We illustrate its use by applying it to the evolution 
of a grain boundary groove in a high vapour pressure material and show that Mullins's classical 
equation describing transport driven by gradients in surface curvature is reproduced asymptotically. 
The microscopic contact angle at the groove root is found to be modified over a thin boundary 
layer, and the apparent contact angle is determined. An explicit transport coefficient is derived that 
governs the evolution rate of systems controlled by surface transport through premelted films. The 
transport coefficient is found to depend on temperature and diverges as the bulk melting temperature 
is approached. 
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It is well established [? ] that melting in any ma- 
terial is initiated at the free surface. Thus, a molecu- 
larly thin layer of melt liquid can exist at the surface at 
temperatures below the bulk transition temperature, de- 
pending on the predominant intermolecular interactions 
at the solid-vapour interface. Typically, repulsive van der 
Waals forces or electrostatic interactions act to thicken 
the film, while the film width is constrained by an attrac- 
tive pressure arising from the liquid being held below its 
freezing point. (For a wide ranging review, see Dash et 
al. Wettlaufer and Worster [2| have shown that it 

is possible to use lubrication theory to model the flow in 
the film, and they found that their results were in good 
agreement with experiment, using bulk parameter values 
(e.g. dynamic viscosity) in the flow |3'J. Furthermore, 
there exists much evidence, both theoretical and 
experimental 0] showing that the continuum fluid dy- 
namical approach, combined with bulk parameters can 
be used for films thicker than about ten molecular diam- 
eters. Therefore, in systems with sufficiently thick films, 
we can derive a system of continuum equations to model 
the evolution of solid-vapour surfaces near to their melt- 
ing point. 

In this letter, we derive the equations pertaining to 
surface melting and flow in a surface-melted film. We 
then apply these to the evolution of a grain-boundary 
(GB) groove and demonstrate that, in the long time limit, 
the original equation derived by Mullins Q 



Vt + By x 



= 0, 



(1) 



is recovered asymptotically. Here, y{x,i) is the surface 
height, x is distance parallel to the surface and B is a con- 
stant transport parameter. We use typical bulk parame- 
ters for the water/ice system at — 1° to give predictions 
of grooving rates that are in agreement with typical ob- 
servations. We focus on ice because it is experimentally 
accessible and because there is much research on it due 
to its unique importance in geoscience. It is also known 
to premelt against air @, and may be doped to yield 



sufficiently thick films to justify the continuum approach 
I- 

Firstly, using the assumption that the system is in local 
thermodynamic equilibrium, we have as a consequence of 
the Gibbs-Duhem relationship [ljj that 



Ps - Pi = Psqn 
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where p, p, T, T m and q m are density, pressure, system 
temperature, melting temperature and latent heat of fu- 
sion respectively, and subscripts s and I correspond to 
the solid and liquid phases. We have also used the as- 
sumption that the system is isothermal. 

Defining h\(x,t) to be the height of the solid-liquid 
surface, h 2 {x, t) to be the height of the liquid-air surface, 
d(x, t) — h>2 — hi to be the thickness of the liquid layer 
(cf Fig. Q]for a specific example), -f s i and -fi a to be the 
(constant) surface free energies of the two interfaces, and 
A to be the effective Hamaker constant [|J , we have that 
the Hclmholtz free energy 



F 



-p s {h 2 -d)- pid + p a h 2 + 7ia(l + /I?) 172 



+ ( 7 ' , + T2^) (1 + ' , ' ) " 2 



dx. 
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Treating the integrand as a function of h\ and d and 
minimising F using the Euler-Lagrange equation gives 



Ps ~Pl 
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where K\ is the curvature of the solid-liquid interface 
z = hi. 

We assume that the thickness of the film is much 
greater than a characteristic molecular dimension, and 
that the slopes of the surfaces are small. Therefore 
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K < 1, 



(5) 
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and we can approximate 



Ps - Pi 
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(6) 



which we combine with Eq. to give a general result 
for a surface-melted film (cf [l|) 



PsQn 



T m — Tn 



= -isiK + 
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(7) 



The thickness of the layer is determined by the com- 
petition between the rise in energy required to sustain 
a liquid film below the melting point, and the van der 
Waals (VdW) forces, as expressed in the balance between 
the LHS and last term on the RHS in ([7]). As can also 
be seen in the equation, the Gibbs-Thomson effect mani- 
fests itself in that an increase in surface curvature causes 
an increase in the thickness of the film. 

The second Euler-Lagrange relationship from above 
yields 



Pi = Pa- llaK - 
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(8) 



So, in regions where the film is thick and has high 
curvature, the liquid pressure is high. This drives flow 
away, towards regions where the curvature is lower and 
the film thinner. 

Finally, assuming that the thickness of the film is suf- 
ficiently small, we can use lubrication theory [11] to give 
the velocity in the layer (see Fig. ([T]) insert) using a no- 
slip condition at the solid-liquid interface and a stress- 
free condition at the liquid-air interface. Mass conserva- 
tion is expressed by 



dq 

Ox 



0. 



(9) 



where q(x) = ~d 3 p x /(3[i), and /x is the dynamic viscosity 
of water. 

We introduce dimensionless variables with a length- 
scale determined by taking the balance between VdW 
and curvature terms in equation ([7]) to give 
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(10) 
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and we write (hi, /12, d, do) = €(%!,%, V, T>q). 

Thus we obtain a coupled pair of dimensionless partial 
differential equations 



-1 - J- - -V" 



V" 



(12) 



n = o, 



(13) 



where 7 = ^yia/jsi- Boundary conditions for these equa- 
tions depend on the specific problem to be addressed. 




FIG. 1: Schematic diagram of a grain-boundary groove. The 
system is isothermal with T — To. 

We can now, for example, address the problem of a 
grain-boundary groove. The model to be analysed is 
shown in Fig. [TJ Two ice crystals, symmetric about the 
GB between them are exposed to the air. The system 
is maintained at a uniform temperature To sufficiently 
close to the melting point T m that surface melting oc- 
curs and a liquid film is formed, separating solid from air. 
Initially, the upper surfaces of both crystals are planar, 
but equilibrium demands a fixed dihedral angle (given by 
the Young-Dupre relationship [Hj]), at the tri-junction 
formed by the GB and the solid-liquid interfaces. The 
premelted film is thereby made thicker near the GB, the 
disjoining force across it is correspondingly weaker and 
the liquid pressure higher. This drives a flow in the pre- 
melted film transporting mass outwards to form a GB 
groove. 

Using this model, we take boundary conditions to 
equations dHJ) , ((T3|) to be D(oo) = £> , "H'(O) = 0, 
T>'(0) = —a and q(0) — 0, corresponding respectively 
to constant film thickness at infinity, constant contact 
angle at the groove base (with cotangent a), continuous 
derivative of the liquid-air interface above the groove and 
zero mass flux across the plane of symmetry. 

As time proceeds, the groove widens and deepens, 
while the thickness of the premelted film stays relatively 
constant. Therefore we can make the assumption (justifi- 
able a posteriori) that T> and hence V" <^ %" away 
from a neighbourhood of the z axis. Thus (fT2j) becomes 
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-H". (14) 
Therefore, assuming that T> w T> (see below) and com- 
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bining with (|13[) . we recover the Mullins equation 

n T +V 3 (j + l)H^ -0, (15) 
with dimensional transport coefficient 

A(rf + 1)7,1 r m 



B 



18irpp s q m (T m - T ) ' 
This equation admits a similarity solution 0] 

H s = [P 3 ( 7 + l)r] 1 / 4 /W 



where 



and the function / satisfies 

f nf 

r+ t_nL = 0t 



(16) 



(17) 



(18) 



(19) 



subject to the boundary conditions f"'(0) = 0, /'(0) = /3, 
/'(oo) = and f(oo) = 0. This yields H s (0) which we 
will use to match to the inner boundary layer solution. 
The constant (3 is the cotangent of the apparent contact 
angle seen by this outer region, and will be determined 
by matching to the inner asymptotic solution, derived as 
follows. 

In the similarity solution regime, we take H ~ T l l 4 f(rj) 
and T> ~ g{rj). Combining these estimates with Eq. (fT2")) 
yields V = Vq + O^t 1 ^ 4 ) and hence the similarity solution 
is expected to break down at the GB when V — > —a. 

In this regime, we assume that the r 1 / 4 time depen- 
dence of the similarity solution is carried over into the 
inner solution. Then consideration of the order of mag- 
nitude of terms in Eq. (fl~3|) shows that H T vanishes for 
large times, and we can make quasi-steady approxima- 
tions of H(£,t) = Ct 1 ' 4 + H(0 and V(£,t) = £>(£) 
inside the boundary layer. 

Applying a zero flux boundary condition at the GB, 
Eq. ([13]) becomes 



D' 



1 H" 



1 

IP 



(20) 



which we integrate with H"(oo) = D"(oo) = 0, D'(0) 
—a and W'(oo) = to obtain 



Hi 



Ct 1 ' 4 + 



D 



c. 
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where, from Eq. (IT21 . D satisfies the equation 



1 



1 



(21) 



(22) 



Substituting H s = r 1 / 4 /(r/) into Eq. ([T3]), we find that 
as t -> oo, TL S satisfies the full equation. Therefore we 



TABLE I: Table of typical values for ice/water/air at T = — 1° 
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1.40 x 10" 26 
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can match Eqs. (IT71) and (|2T|) by letting Hi nne r{0 ^ H s 
and -D(£) — > T>q for large £ (the outer edge of the inner 
solution), and matching this to the small £ limit of Eq. 

(US}. 

For small £, Eq. (fl~7| becomes 

U s =H s (0)+£/? + O(t- 1 / 4 ), 



giving that 



Hi 



D-V 



(23) 



(24) 



1+7 1+7 
Importantly, we have a value for the apparent cotangent 

P = T^~- (25) 



1 
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Eq. (I2"5li shows that for 7; a large relative to 7 S ; , the sys- 
tem will try to minimize the liquid-air interfacial area, 
and hence /3 — > 0. Conversely for 7 S ; 3> 7; OJ /3 — > a 
as the solid-liquid interfacial curvature is minimized at 
the GB to reduce surface area. This also implies that 
the apparent cotangent /3 can differ substantially from 
the actual cotangent a suggesting that experimental de- 
termination of a could yield incorrect results due to the 
apparent contact angle being measured instead. 

To solve the complete system of equations, we used an 
explicit numerical scheme. At each step, the T> profile is 
calculated from (fT2"]) by a relaxation technique. The new 
value of H at (r + dr) can then be calculated from (fl3|) . 
using the previous values of %{t) and T>(t). 

Fig. [2] shows the time evolution of the groove. Good 
agreement can be seen, even at small times, between the 
numerical solutions (continuous lines) and the similarity 
solution (dashed lines). Boundary conditions for the sim- 
ilarity solution are H'(0) = (3, q(0) — 0, H(oo) = and 
"W(oo) = 0. The inset figure shows the inner boundary 
layer region for t — 2.9 x 10~ 4 . The inner asymptotics 
(dash-dot) are shown, and also show good agreement 
with the numerical solution. 
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FIG. 2: Numerical solutions (continuous lines) and similarity 
solutions (dashed lines) for time evolution of a groove at a 
liquid-air interface, y is the the height of the interfaces rela- 
tive to the initial position. Graphs are plotted in terms of x 
and t at t = 8.7 x 10~ 7 , 1.7 x 10" 5 , 2.9 x 10" 4 . Inset shows nu- 
merical solution (continuous line), similarity solution (dashed 
line) and inner boundary layer asymptotics (dash-dot line) for 
the inner boundary layer regime at t = 2.9 x 10 -4 . Parameters 
were chosen as "Do = 1, 7 = 0.5. 

From Eq.(fl~6f we see that the rate-controlling coeffi- 
cient is inversely proportional to the reduced tempera- 
ture, and grows as To approaches the melting point. We 
also see that less viscous, less dense liquids with thicker 
layers (A large) and stronger surface tension cause faster 
grooving than more viscous, dense liquids with thin lay- 
ers and weak surface tension. Typical parameters for 
ice at -1° (see Tabic [TJ) give B to be approximately 
1.4 x 10 -26 m 4 s -1 and so typical lengthscales (propor- 
tional to (Bt 1 / 4 )) are of order 0.35/im in Is and 3.5/um 
in 3 hours, which is in line with observations of ice [13j . 
At this time, the authors are unaware of any experimen- 
tal measurements of B in any material within a suitable 
temperature range. 

The analysis presented above is relevant to any mate- 
rial that premelts against its atmosphere and forms a con- 
stant contact angle at the groove root, in circumstances 
where it is maintained at a suitable temperature that 
premelting can occur. However, several points should be 
noted, (i) We have assumed that the slope of the inter- 
faces is always small, and therefore some detail of the 
form of the groove root may have been lost, (ii) We as- 
sumed that there is a complete phase change across the 
surface z = Hi, while the transition may actually occur 
across several atomic diameters and involve some level 
of ordering(eg [? ]). However, provided the premelted 
film is thick enough, the affects of this should be negli- 
gible and we may use continuum approximations. Ap- 
propriately thick films may be ensured by the addition 
of dopants to the system [§]. (iii) We have considered 
surface melting controlled by unretarded VdW forces. 
This can be easily generalised to retarded VdW forces 



and electrostatic interactions [3j. (iv) We require a to 
be small for the lubrication approximation to be valid in 
the inner region. However, groove root angles for ice may 
be very small, so this analysis can break down near the 
groove root. 

So, in conclusion, we have demonstrated a new mech- 
anism for mass transport on the surface of a crystal in 
a suitably high temperature regime. Surface melting is 
a ubiquitous phenomena, and so the equations that we 
have derived should be applicable to a wide range of ma- 
terials and situations. We have illustrated the mecha- 
nism of surface transport by applying it to the case of a 
grain-boundary groove and demonstrated that it yields 
Mullins's equation in the asymptotic limit, although the 
dihedral angle determined by the surface energies is mod- 
ified over the boundary layer by a factor 1/(1 +7). We 
note that much richer behaviour will result from effects 
such as crystal anisotropy, impurity content, asymmetry 
and other grooving mechanisms Q ■ 

The most important result of this letter is that we have 
found an explicit transport parameter B (Eq. (|16p) that 
controls the rate of surface transport in the system. It 
should be noted that B docs not depend on the geometry 
of the grain-boundary groove. Therefore this parameter 
should be relevant to most problems involving surface 
melting such as sintering, grain-boundary groove migra- 
tion, crystal growth and annealing to name but a few. 

The authors would like to thank J.S.Wettlaufer for his 
critical readings of the letter. 
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